R Notebook

只分析ctrl

source("tianfengRwrappers.R")

Ctrl

umap

Return:

umap_plot

markers overview

Return:

markers_heatmap

ctrl.markers <- FindAllMarkers(ctrl, min.pct = 0.1, 
                                    logfc.threshold = 0.5, only.pos = T)
top5markers <- ctrl.markers %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_logFC)
markers_heatmap <- dhm(top5markers$gene, ctrl)
(markers_heatmap + theme(axis.text = element_text(size= 14, color = "black")))%>%
  ggsave(filename = "./ctrl_only/ctrl_markergene.png", device = png, height = 9, width = 16, plot = .) 

dotplot

top3markers <- ctrl.markers %>% group_by(cluster) %>% slice_max(n = 3, order_by = avg_logFC)
# dotplot <- DotPlot(ctrl, features = top3markers$gene, cols = c("lightgrey","#ff2121")) + coord_flip() + RotatedAxis() + theme_classic() #用来获得白色 而不是透明的背景

dotplot <- Dotplot(top3markers$gene[!duplicated(top3markers$gene)], ctrl) + coord_flip()

##单个基因featureplot ### Return: gene_plot

gene_plot <- f("Gasp",ctrl)
#ggsave(filename = paste0(gene_name,".svg"), device = svg, height = 5, width = 6, plot = gene_plot) #对svg调整坐标轴

##两个基因umap共定位 ### Return: doublegene_featureplot

gene_name <- c("Tom","peb")
p <- doublegene_featureplot(gene_name, ctrl)
ggsave(filename = paste0(gene_name[1],"vs",gene_name[2],".png"), device = png, height = 5, width = 20, plot = p) #对svg调整坐标轴

多基因堆叠小提琴图

Return:

violin_plot

violin_plot(c("salm","kni","Dl","ct","wg","vn","bs"),ctrl) %>%
  ggsave(filename = paste0("violinplot.png"), device = png, height = 5, width = 5, plot = .)

多图featureplot

Return:

mulit_featureplot

gene_name <- c("bnl", "btl", "htl", "Egfr") #可视化的基因名称
multi_featureplot(gene_name, ctrl) %>% ggsave(paste0("./ctrl_only/",as.character(gene_name), ".png"),device = png, width = 6, height = 5, plot = .)

multi_featureplot(c("upd2","dome","PlexA","PlexB"), ctrl) %>% ggsave("./ctrl_only/upd2dome.png",device = png,height = 6, width = 7, plot = .)

multi_featureplot(c("Sema5c","Sema1a","Sema2a","Sema1b"), ctrl)%>% ggsave("./ctrl_only/sema.png",device = png,height = 6, width = 7, plot = .)

multi_featureplot(c("peb","Dad","dpp","ex"),ctrl) %>% ggsave("./ctrl_only/pebDadex.png",device = png,height = 6, width = 7, plot = .)

转录因子 SCENIC

Return:

TFheatmap

regulonActivity <- read.csv("./regulonActivity.csv",row.names = 1)
topTFs <- read.csv("./topRegulators.csv")
top5 <- topTFs %>% group_by(seuratCluster) %>% slice_max(n = 5, order_by = RelativeActivity)

#根据已知的top5更新行顺序
regulonActivity <- regulonActivity[sapply(top5$Regulon, function(e) {which(rownames(regulonActivity) == e)}), ]
#annotation_col <-  data.frame(cluster = factor(top5$seuratCluster), row.names = make.names(top5$Regulon,TRUE)) #make.names用来生成不冲突的行
regulonActivity <- regulonActivity[,c("DB","DT","NE","PC","SB","TC","VB")]

TFheatmap <- pheatmap(regulonActivity, breaks = unique(c(seq(-3,3, length=400))), 
                 color = colorRampPalette(c("#1E90FF", "white", "#ff2121"))(400),
                border_color = NA, cluster_rows = FALSE, cluster_cols = FALSE,
                main = "regulonActivity",angle_col = 45, show_rownames = T)
#ggsave(filename = paste0("regulonActivity",".svg"), device = svg, height = 8, width = 12, plot = TFheatmap)

pseudotime featureplot

Return:

pseudotime in featurePlot

pseudotimePlot <- f("pseudotime", cols = c("hotpink","#1E90FF"))
#ggsave(filename = paste0("pt_featureplot2",".svg"), device = svg, height = 5, width = 6, plot = pseudotimePlot)

基因+拟时

Return:

genewithpseudotime

datMat <- as.matrix(GetAssayData(ctrl,slot = "data"))
df <- data.frame(t(ctrl$pseudotime))
df <- do.call(data.frame,lapply(df, function(x){
  x = x/20
  replace(x, is.infinite(x),0)}))
rownames(df) <- "pseudotime"
colnames(df) <- colnames(datMat)
datMat <- rbind(df, datMat)

test <- CreateSeuratObject(counts = datMat)
test@reductions[["umap"]] <- ctrl@reductions[["umap"]]
test@active.ident <- ctrl@active.ident
# SetAssayData(ctrl,slot = "data", new.data = datMat)
#共定位
gene_name <- c("Tom","pseudotime")
genewithpseudotime <- FeaturePlot(test, pt.size = 0.8, slot = "counts",min.cutoff = 0, max.cutoff = 4, features = gene_name, cols = c("lightgrey","#1E90FF", "#ff2121"), label = TRUE, blend = T, blend.threshold = 0) + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())
#ggsave(filename = paste0(gene_name[1],"vs",gene_name[2],".png"), device = png, height = 5, width = 20, plot = genewithpseudotime)

基因拟时表达

Return:

pseudotime_cells pseudotime_bins pseudotime_ridges pseudotime_ridges2

#选取细胞亚群进行分析
cell.subset <- subset(x = ctrl,idents = c("DT","SB"),invert = FALSE )
nbins <-  15
cell.subset <- make_bins(cell.subset,nbins)
top_specific_marker_ids <- readRDS("top_specific_marker_ids.RDS")
pseudotime_cells <- make_pseudotime_heatmap_bycells(cell.subset,top_specific_marker_ids,50)
#ggsave(filename = "pseudotime_heatmap_bycells.png", device = png, height = 6, width = 10, plot = pseudotime_cells)

pseudotime_bins <- make_pseudotime_heatmap_average(cell.subset,top_specific_marker_ids,50)
#ggsave(filename = "pseudotime_heatmap_bins.png", device = png, height = 6, width = 10, plot = pseudotime_bins)

df <- get_draw(ctrl, genelist = head(top_specific_marker_ids,10), expr_threshold = 1)

pseudotime_ridges <- ggplot(df, aes(x = pseudotime , y = gene , fill = gene)) +
 ggridges::geom_density_ridges(alpha = 0.5,show.legend =T) + labs(title = 'DEGs with pseudotime')+
 theme_ipsum() +  theme(legend.position="none", panel.spacing = unit(0.1, "lines"),
                        strip.text.x = element_text(size = 8))
#ggsave(filename = "pseudotime_ridges.png", device = png, height = 8, width = 10, plot = pseudotime_ridges)

pseudotime_ridges2 <- ggplot(df, aes(x = pseudotime, y = gene, fill = ..x..))  +
  geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01,alpha = 0.5) +
  scale_fill_viridis(name = ctrl$pseudotime, option = "C") +
  labs(title = 'DEGs with pseudotime') + theme_ipsum() +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      strip.text.x = element_text(size = 8))
##ggsave(filename = "pseudotime_ridges2.jpg", device = jpeg, height = 8, width = 10, plot = pseudotime_ridges2)

GO

Return:

GO_plot

library(org.Dm.eg.db)
for (cl in levels(Idents(ctrl))){
  #选择cluster对应的marker
  gene_list <- ctrl.markers[ctrl.markers$cluster==cl,]$gene
  #ID转换
  genes <- bitr(gene_list, fromType="SYMBOL",toType=c("ENTREZID","ENSEMBL"),
                   OrgDb = org.Dm.eg.db)
  enrich.go <- enrichGO(gene = genes$ENTREZID, #基因列表文件中的基因名称
                          OrgDb = org.Dm.eg.db,
                          keyType = 'ENTREZID',
                          ont = 'ALL', #可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
                          pAdjustMethod = 'fdr',
                          pvalueCutoff = 0.05, 
                          qvalueCutoff = 0.2, 
                          readable = TRUE)
  
  GO_plot <- dotplot(enrich.go,title = paste(cl,"GO"),showCategory = 15) + theme_classic() + theme(text = element_text(colour = 'black', size=12))
  ggsave(filename = paste0("./ctrl_only/GO",cl,"_GO.svg"), device = svg, height = 8, width = 10, plot = GO_plot)
}
# barplot(enrich.go,showCategory = 15) + theme_classic()
# cnetplot(enrich.go)

cellchat

Return:

act_weight chord_plot1 chord_plot2 pathway_heatmap pathway_chord act_bubble1 act_bubble2

cellchatobj <- readRDS("cellchatobj.rds")
groupSize <- as.numeric(table(cellchatobj@idents))
# par(mfrow = c(1,2), xpd=TRUE)
#netVisual_circle(cellchatobj@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
act_weight <- netVisual_circle(cellchatobj@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength",margin = 0.2)

#chord plot
# par(mfrow = c(1,2), xpd=TRUE)
chord_plot1 <- netVisual_chord_gene(cellchatobj, sources.use = c(6), targets.use = c(1:7), lab.cex = 1, legend.pos.y = 30, legend.pos.x = 10)
chord_plot2 <- netVisual_chord_gene(cellchatobj, sources.use = c(1:7), targets.use = c(6), lab.cex = 1,legend.pos.y = 30, legend.pos.x = 10)


pathways.show <- c("HEDGEHOG signaling pathway") 
# Hierarchy plot

vertex.receiver = seq(1,5) # a numeric vector. 
# netVisual_aggregate(cellchatobj, signaling = pathways.show, layout = "hierarchy",vertex.receiver = vertex.receiver)
# netVisual_aggregate(cellchatobj, signaling = pathways.show, layout = "circle")
pathway_chord <- netVisual_aggregate(cellchatobj, signaling = pathways.show, layout = "chord")


pathways.show <- c("NOTCH signaling pathway") 
par(mfrow=c(1,1))
pathway_heatmap <- netVisual_heatmap(cellchatobj, signaling = pathways.show, color.heatmap = c("#1E90FF","#ff2121")) 
#cellchatobj@netP 对这里的元素画图
pathway_heatmap

for(i in c(1:7)){
act_bubble1 <- netVisual_bubble(cellchatobj, sources.use = i, targets.use = c(1:7), remove.isolate = FALSE, grid.on = TRUE)+ theme_classic() + geom_point(size = 2) +
  theme(text = element_text(family = NULL, colour = 'black', size=10)) #6 PC作为source
act_bubble2 <- netVisual_bubble(cellchatobj, sources.use = c(1:7), targets.use = i, remove.isolate = FALSE, grid.on = TRUE) + theme_classic() + geom_point(size = 2) + 
  theme(text = element_text(family = NULL, colour = 'black', size=10)) #PC作为target
png(filename = paste0(cellchatobj@idents[i],".png"),width = 3500, height = 1800, res = 500)
print(act_bubble1)
dev.off()
png(filename = paste0(cellchatobj@idents[i],"2.png"),width = 3500, height = 1800, res = 500)
print(act_bubble2)
dev.off()
}

组合图片

fig.1left <- plot_grid(umap_plot, markers_heatmap, labels = "AUTO", ncol = 1, nrow = 2, 
                   label_size = 18, rel_heights = c(1, 1), rel_widths = c(1, 1), align = "v")
fig.1 <- plot_grid(fig.1left, dotplot, labels = c("","C"), ncol = 2, nrow = 1, 
                   label_size = 18, rel_heights = c(1, 1),rel_widths = c(1.1, 0.8), align = "v") + theme_void()
ggsave(filename = paste0("fig1r",".png"), device = png, height = 10, width = 16, plot = fig.1)



fig.2 <- mulit_featureplot
ggsave(filename = paste0("fig2",".png"), device = png, height = 10, width = 12, plot = fig.2)


fig.3 <- plot_grid(act_weight, pathway_chord, act_bubble1, act_bubble2, labels = "AUTO", ncol = 2, nrow =2, 
                   label_size = 18,align = "hv" )
ggsave(filename = paste0("fig3",".png"), device = png, height = 10, width = 16, plot = fig.3)

fig.4 <- plot_grid(corheatmap$gtable, labels = "AUTO", label_size = 18, align = "hv" ) + theme_classic() + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())
ggsave(filename = paste0("fig4",".png"), device = png, height = 8, width = 8, plot = fig.4)


fig.5upper <-  plot_grid(subumap_plot, DBsubumap, labels = c("A","B"), label_size = 18, ncol = 2, nrow = 1)
fig.5 <- plot_grid(fig.5upper, DB_features, labels = c("",""), rel_widths = c(1, 1), label_size = 18, ncol = 1, nrow = 2)
ggsave(filename = paste0("fig5",".png"), device = png, height = 10, width = 10, plot = fig.5)

fig.6 <-  plot_grid(pseudotime_bins$gtable, pseudotime_ridges, pseudotime_ridges2, labels = "AUTO", ncol = 1, nrow = 3, 
                   label_size = 18) + theme_classic() + theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())
ggsave(filename = paste0("fig6",".png"), device = png, height = 18, width = 10, plot = fig.6)

HSD

umap

umap_plot <- umapplot(HSD, group.by = "celltype")

Scale for 'colour' is already present. Adding another scale for 'colour', which will
replace the existing scale.

umap_plot
ggsave(filename = "./HSD_only/HSD_umap.png", 
       device = png, height = 6, width = 7, plot = umap_plot)

markers overview

Return:

markers_heatmap

HSD.markers <- FindAllMarkers(HSD, min.pct = 0.1, 
                                    logfc.threshold = 0.5, only.pos = T)
top5markers <- HSD.markers %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_logFC)
markers_heatmap <- dhm(top5markers$gene, HSD)
(markers_heatmap + theme(axis.text = element_text(size= 14, color = "black")))%>%
  ggsave(filename = "./HSD_only/HSD_markergene.png", device = png, height = 9, width = 16, plot = .) 

dotplot

top3markers <- HSD.markers %>% group_by(cluster) %>% slice_max(n = 3, order_by = avg_logFC)
## dotplot <- DotPlot(HSD, features = top3markers$gene, cols = c("lightgrey","##ff2121")) + coord_flip() + RotatedAxis() + theme_classic() ##用来获得白色 而不是透明的背景

dotplot <- Dotplot(top3markers$gene[!duplicated(top3markers$gene)], HSD) + coord_flip()

##两个基因umap共定位 #### Return: doublegene_featureplot

gene_name <- c("Tom","peb")
p <- doublegene_featureplot(gene_name, HSD)
ggsave(filename = paste0('./HSD_only/',gene_name[1],"vs",gene_name[2],".png"), device = png, height = 5, width = 20, plot = p)

多基因堆叠小提琴图

Return:

violin_plot

violin_plot(c("salm","kni","Dl","ct","wg","vn","bs"),HSD) %>%
  ggsave(filename = paste0("./HSD_only/HSD_violinplot.png"), device = png, height = 5, width = 5, plot = .)

多图featureplot

Return:

mulit_featureplot

gene_name <- c("bnl", "btl", "htl", "Egfr") ##可视化的基因名称
multi_featureplot(gene_name, HSD) %>% ggsave(paste0("./HSD_only/",as.character(gene_name), ".png"),device = png, width = 6, height = 5, plot = .)

multi_featureplot(c("upd2","dome","PlexA","PlexB"), HSD) %>% ggsave("./HSD_only/upd2dome.png",device = png,height = 6, width = 7, plot = .)

multi_featureplot(c("Sema5c","Sema1a","Sema2a","Sema1b"), HSD)%>% ggsave("./HSD_only/sema.png",device = png,height = 6, width = 7, plot = .)

multi_featureplot(c("peb","Dad","dpp","ex"),HSD) %>% ggsave("./HSD_only/pebDadex.png",device = png,height = 6, width = 7, plot = .)

转录因子 SCENIC

Return:

TFheatmap

regulonActivity <- read.csv("./regulonActivity.csv",row.names = 1)
topTFs <- read.csv("./topRegulators.csv")
top5 <- topTFs %>% group_by(seuratCluster) %>% slice_max(n = 5, order_by = RelativeActivity)

##根据已知的top5更新行顺序
regulonActivity <- regulonActivity[sapply(top5$Regulon, function(e) {which(rownames(regulonActivity) == e)}), ]
##annotation_col <-  data.frame(cluster = factor(top5$seuratCluster), row.names = make.names(top5$Regulon,TRUE)) ##make.names用来生成不冲突的行
regulonActivity <- regulonActivity[,c("DB","DT","NE","PC","SB","TC","VB")]

TFheatmap <- pheatmap(regulonActivity, breaks = unique(c(seq(-3,3, length=400))), 
                 color = colorRampPalette(c("##1E90FF", "white", "##ff2121"))(400),
                border_color = NA, cluster_rows = FALSE, cluster_cols = FALSE,
                main = "regulonActivity",angle_col = 45, show_rownames = T)
##ggsave(filename = paste0("regulonActivity",".svg"), device = svg, height = 8, width = 12, plot = TFheatmap)

GO

Return:

GO_plot

library(org.Dm.eg.db)
for (cl in levels(Idents(HSD))){
  ##选择cluster对应的marker
  gene_list <- HSD.markers[HSD.markers$cluster==cl,]$gene
  ##ID转换
  genes <- bitr(gene_list, fromType="SYMBOL",toType=c("ENTREZID","ENSEMBL"),
                   OrgDb = org.Dm.eg.db)
  enrich.go <- enrichGO(gene = genes$ENTREZID, ##基因列表文件中的基因名称
                          OrgDb = org.Dm.eg.db,
                          keyType = 'ENTREZID',
                          ont = 'ALL', ##可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
                          pAdjustMethod = 'fdr',
                          pvalueCutoff = 0.05, 
                          qvalueCutoff = 0.2, 
                          readable = TRUE)
  
  GO_plot <- dotplot(enrich.go,title = paste(cl,"GO"),showCategory = 15) + theme_classic() + theme(text = element_text(colour = 'black', size=12))
  ggsave(filename = paste0(cl,"_GO.svg"), device = svg, height = 8, width = 10, plot = GO_plot)
}
## barplot(enrich.go,showCategory = 15) + theme_classic()
## cnetplot(enrich.go)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.